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Abstract 

In the 1960s, a moment approach to linear time varying (LTV) minimal norm 
impulsive optimal control was developed, as an alternative to direct approaches 
(based on discretization of the equations of motion and linear programming) or in- 
direct approaches (based on Pontryagin's maximum principle). This paper revisits 
these classical results in the light of recent advances in convex optimization, in par- 
ticular the use of measures jointly with hierarchy of linear matrix inequality (LMI) 
relaxations. Linearity of the dynamics allows us to integrate system trajectories 
and to come up with a simplified LMI hierarchy where the only unknowns are mo- 
ments of a vector of control measures of time. In particular, occupation measures 
of state and control variables do not appear in this formulation. This is in stark 
contrast with LMI relaxations arising usually in polynomial optimal control, where 
size grows quickly as a function of the relaxation order. Jointly with the use of 
Chebyshev polynomials (as a numerically more stable polynomial basis), this allows 
LMI relaxations of high order (up to a few hundreds) to be solved numerically. 



1 Introduction 

In the 1960s, it was realized that many physically relevant problems of optimal control 
were inappropriately formulated in the sense that the optimum control law (a function 
of time and/or state) cannot be found if the admissible functional space is too small. 
This observation was the main driving force of the papers (T7J [18], [20] which introduce 
optimal control problems formulated in the space of measures. The approach is well 
summarized in the textbook [15], which contains some (academic) examples of optimal 
control problems without solutions. This motivated the introduction of many concepts 
of functional analysis (density, completeness, duality, separability) in control engineering, 
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building up on the advances on mathematical control theory and calculus of variations. 
As promoted in [15], an optimal control problem should be formulated in the dual of a 
Banach space which is large enough for the solution to be attained. 

Most of the literature on numerical optimal control focuses on direct approaches (based on 
discretization of the equations of motions and linear programming) or indirect approaches 
(based on the necessary optimal conditions of Pontryagin's maximum principle). When 
applied to optimal control problems whose optima cannot be attained, these numerical 
approaches typically face difficulties. In this context, we believe that it is timely to revisit 
classical results by Neustadt [17] on the formulation of optimal control problems for linear 
time varying (LTV) systems as a problem of moments, where the decision variables (from 
which an optimal control law can be extracted) are measures subject to a finite number 
of linear constraints. 

There is an important literature, especially from the 1960s, on moment formulations to 
optimal control of ordinary differential equations (ODEs) and partial differential equa- 
tions (PDEs), see e.g. [3] and references therein, as well as the comments of [?, p. 586]. 
Currently, this approach is not frequently used by engineers, and in our opinion this may 
be due, on the one hand, to the technicality of the underlying concepts of functional 
analysis, and, on the other hand, to the absence of numerical methods to deal satisfacto- 
rily with optimization problems in large functional spaces such as spaces of measures or 
distributions. Regarding the first point, we strongly recommend the textbook [15] which 
is a very readable account of elementary functional analysis useful for engineers. Regard- 
ing the second point, there has been recent advances in convex optimization, especially 
semidefinite programming (optimization over linear matrix inequalities, LMIs), for solving 
numerically generalized problems of moments, i.e. linear programming (LP) problems in 
Banach spaces of measures, see p2J E] and references therein. 

Our contribution is therefore to revisit the classical formulation by Neustadt [T7] in the 
light of recent advances on LMI hierarchies for solving generalized problems of moments. 
In our previous work [11] , we formulated polynomial optimal control problems with semi- 
algebraic state and control constraints as generalized problems of moments that can be 
solved with asymptotically converging LMI hierarchies. The optimal control problem 
is relaxed to an LP problem in the space of occupation measures, which are measures 
of time, state and control encoding the trajectories of the system. The main drawback 
of this approach is the rapid growth of the size of the LMI problems in the hierarchy, 
making the approach applicable to small-size problems only (say at most 3 states and 
2 controls). In the current paper, which focuses on the specific case of LTV dynamics 
depending affinely on the control variable, we first replace control variables by interpret- 
ing them as measures of time. This is similar to our previous work [6] which also dealt 
with impulsive control design in a more general setting. Second, we get rid of the state 
variables by integrating numerically the LTV ODE. This is possible because the ODE 
depends linearly on the state and the control. As a result, the optimal control problem 
is relaxed to an LP on measures depending only on time. It follows that there is no need 
for an LMI hierarchy since in the univariate case finite-dimensional moment LMI condi- 
tions are necessary and sufficient. However, there is still a hierarchy of LMI conditions 
in connection with the polynomial approximations of increasing degree we use to model 
the integrated system trajectories. To deal with those high degree univariate polynomials 
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and moment matrices, we use Chebyshev polynomials instead of monomials. Indeed, high 
degree Chebyshev polynomials behave much better numerically than monomials, and we 
can rely on functionalities of the chebfun package for Matlab [221123] to integrate the LTV 
ODE and manipulate polynomials. To illustrate the above methodology, we show on some 
examples how to obtain very good approximations of impulse times and amplitudes of an 
optimal solution. 



2 Relaxed linear optimal control 

Consider the linear time varying (LTV) optimal control problem 

"l r-tp 

q* := inf \\u\\i := / \v,j(t)\dt 

s.t. x(t) = A(t)x(t) + B{t)u(t) (!) 
x(ti) G lR n given 
x(£f) £ R n given 

where the minimization is w.r.t. a vector of control functions Uj G Z^Qij, tp\), j = 
1, . . . , m, and A G £°°([tj, t F ]; R nxn ), B G <7([tj, R nxm ), on a given bounded time 
interval [tj,tp} C R. 

In general the infimum is not attained, and the optimal control problem is relaxed to 

p* := inf \\fi\\ TV :=^2 \Hj\(di) 

i=i Jtl 

s.t. x(dt) = A(t)x(t)dt + B(t)/i(dt) ( 2 ) 
x(ti) G lR n given 
x(tp) G lR n given 

where the minimization is w.r.t. a vector of (signed) measures fij G M([ti,tp]), j = 
1, . . . , m, and ||/x||tv denotes the total variation, or norm, of vector measure /i. A measure 
in M([ti, tj?]) of finite norm is identified (by a representation theorem of F. Riesz, see e.g. 
[T9| Section 21.5]) as a continuous linear functional acting on the space of continuous 
functions C([t/,tp]). 

Problem Q is a relaxation of problem ([T]) since we enlarge the space of admissible con- 
trols. Indeed, problem (II]) is equivalent to problem ^ restricted to measures which 
are absolutely continuous w.r.t. time, i.e. fij(dt) = Uj(t)dt for some Uj G L l ([ti, tp]), 
j — 1, . . . , m. The motivation for introducing relaxed problem ^ is as follows. 

Lemma 1 The infimum is attained in problem (flU and it is equal to the infimum of 
problem |l]), i.e. q* = p* . 

The proof of Lemma [T] is relegated to the end of Section [3j We will introduce a numerical 
method to deal directly with relaxed problem ^ in measure space M, bypassing the 
potential difficulties coming from the fact that the infimum in problem is typically 
not attained in function space L l . Before this, we need to reformulate optimal control 
problem ([21) problem of moments. 
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3 Problem of moments 



Now we integrate the differential equation in problem ^ to obtain an equivalent problem 
of moments. For more details, see e.g. [51 Section 2.2]. 

Let fi G W 1 ' 1 ([tj,tF];M. n ) denote the absolutely continuous solution of the Cauchy prob- 
lem fi(t) = A(t)fi(t) with fiiti) equal to the i-th column of I n , the n-by-n identity matrix, 
for i = 1, . . . , n. The matrix 

F(t) := [ h(t) ■■■ f n (t) ] G W^dtjM^n 

therefore satisfies^ the matrix ODE 

F(t) = A(t)F(t), F(ti) = I n , te [tj,t F ]. 

From [5J Theorem 2.2.3j^J matrix is differentiable and any function x[t) satisfying 

x{t F ) = F{t F ) x(ti)+ / F- l {t)B{t)^{dt) 

is a solution to the differential equation 

x(dt) = A{t)x{t)dt + B(t)fi(dt), t G [tj, t F \. (3) 

Letting 

G{t) := (F-\t)B(t)f 

= [ 9x(t) ■■■ 9n(t) ] eC([tj,t F ];R mxn ) } 

h := F~ 1 (t F )x(t F ) - z(tj) G R m , 
we can replace the differential equation ^ with the integral equation: 

f F G T (t)fi(dt) = h. 
Jti 

It follows that problem (|2]) can be written equivalently as 

(4) 



p = mm \\fi\\TV 

s.t. (g h /i) = hi, i = l,...,n 

where 

i*tp m i*tp 

Jti _j Jtj 

denotes the duality bracket between C(\pi,t F ]; M. m ) and its dual M([ti,t F ); IR m ), a bilinear 
form pairing C([t/,t F ];M m ) and M([t h t F ]; R m ). Problem Q is a problem of moments 
consisting of finding m measures subject to n linear constraints. 



: A function belongs to the Sobolev space W 1 ' 1 of absolutely continuous functions if it is the integral 
of a function of Lebesgue space L l , 

2 In this reference the authors define the fundamental matrix solution as a bivariate matrix M(t, s) such 
that dM(t,s)/dt = A(t)M(t,s), M(s,s) = I n . The connection with our definition is that F[t) = M(t,ti). 
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Proof: (of Lemma [Ip The proof that the infimum of problem (|2J) is attained follows from 
[TT1 Theorem 1]. Moreover, in [TTJ Theorem 4] it is shown that there exists a solution 
IX to problem Q, hence to problem ([2]), which is the (signed) sum of at most n Dirac 
measures. Finally, it is shown in [T7J pp. 45-46] that there is a sequence of functions u k G 
ti?]; M m ), k = 1,2,... which are admissible for problem Q, i.e. fi(dt) = u k (t)dt 
satisfies the constraints in problem Q, and which are such that lim^oo = p*. □ 



4 Primal and dual conic LP 

By decomposing each signed measure fij as a difference of two nonnegative measures 
(using the Jordan decomposition theorem, see e.g. [TjJl Section 17.2]), i.e. 



m, 



problem Q can be written as a linear programming (LP) problem on the cone of non- 
negative measures 

p* = inf (l,/i + > + (l,/0 

s.t. (g h - (g h fi ) = hi, i = l,...,n (5) 
/i + > 0, > 

where 1 denotes the m- dimensional vector of functions identically equal to one, and the 
above minimization is w.r.t. two vector- valued nonnegative measures fi + G M([ti, tp]] M m ), 
\i~ G M([tj, tp}] M m ). It is easy to show that problem ^ is equivalent to problem Q. 

Problem ^ is the dual of the following LP on the cone of nonnegative continuous functions 
(see [21] for details of the derivation): 



d* = sup y^^jhj 

i=i 

n 

s.t. z+(t) -=l + Yl Vi9i(t) > * e [*J. ^1 (6) 

n 

z-(t) :=1-J2yi9i(t) >0 te[t/,t F ] 



»=i 



where the maximization is w.r.t. a vector y G M. n parametrizing two vector-valued non- 
negative continuous functions z + G C([tf, t F ]; M m ), z~ G C([tf, M m ). Denoting 

Halloo := sup \zj{t)\ 

tE[ti,tp] ,j=l,...,m 

for any z G C([t 7 ,t F ];M m ), remark that LP problem @ can be also written as 

d* = sup /i T y , , 

s.t. ||G(t)y||oo < I- l J 

Lemma 2 There is no duality gap between LP |5p and |^) ; ie. p* = d* . 



5 



Proof: Define the vector r(fi + ,fi ) G IR n+1 with entries ro(/i + ,/i ) := + ), 

ri(/i + , /i") := - /i - ), i = l,...,n and the set R := {r(//+, /i") : /i") G 

Mi m } C M n+1 where M+ m denotes the cone of nonnegative measures in M([tj, tp); M. 2m ). 
Let us invoke [H Theorem 3.10] which states that p* = d* provided p* is finite and R is 
closed. Finiteness of p* follows immediately since we minimize the total variation. To 
prove closedness, we have to show that all accumulation points of any sequence r(/i+, 
belong to R. Since the supports of the measures are compact and p* is finite, hence 
/i + and /i~ are bounded, the sequence (/z+,/i~) is bounded. By the weak-* compactness 
of the unit ball in the Banach space of bounded signed measures with compact support 
(Alaoglu's Theorem, see e.g. [TSl Section 5.10] or [191 Section 15.1]), there is a subsequence 
(/i+ ft ,/i~ fe ) that converges weakly-* to an element G M^ m . As 1 and all gi belong 

to C([tj,t F ];R m ) then lim^ r(/x+ , //" ) = r(/i + ,/i") G R. □ 

Zero duality gap implies that any optimal pair solving LPs (J5]|6]) sat- 

isfies the complementarity conditions 

(z+, (if) = 0, (Zj~,Vj) = 0, j = 1, . . . , m. 

This means that the support of each measure resp. is included in the set {£ G 
[tj, t F ) : Z +(t) = 0}, resp. {t G [t 7 , t F ) : z"(t) = 0}, for j = 1, . . . , m. 

Note that formulation of the dual problem dates back to the work of Neustadt [T7] 
and was preferred to the primal formulation for numerical solution of the optimal control 
problem. One objective of this paper is to show that the recent advances on the moment 
problem give an efficient computation procedure for the primal problem. In particular, 
one obtains very good approximations of an optimal solution of ^ (impulse times and 
amplitudes). 

Remark 1 The vector G(t)y involved in LP problem |6j) is known as the primer vector 
introduced in the seminal work [13]. This primer vector is defined as the velocity adjoint 
vector arising by the application of Pontryagin's maximum principle to optimal trajectory 
problems. The primer vector must satisfy Lawden's well-known necessary conditions for 
an optimal impulsive trajectory. 



5 Integration of the LTV ODE 

In order to compute matrix F(t), the last step before obtaining a tractable problem, we 
have to integrate numerically the ODE x(t) = A(t)x(t). Matrix G(t) solves numerically 
the LTV system of equations F(t)G T (t) = B{t). In practice, we use the Matlab software 
package Chebfun [22] to build Chebyshev polynomial interpolants of the problem data 
A(t) and B(t). Some of its specialized numerical routines to compute F(t) and G(t) have 
been used here. 

We want to characterize the error of approximating G(t) by its Chebyshev polynomial 
interpolant on d points, see [23] for an introduction on this subject. Define Gd{t) as the 
Chebyshev interpolant of G(t) on d points. Assume furthermore that the error induced 
by algebraic manipulations for obtaining Gd{t) can be properly controlled. Then the 
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approximation error e d := \\G d — G||oo for large d is conditio nned by the regularity of 
G(t). We recall the main results of approximation theory which can all be found in [23] : 

• if G(t) is fc-times continuously differentiable, then e d — > at the algebraic rate of 
0(d~ k ) as d ->• oo; 

• if is fc-times differentiable with its k-th derivative of bounded variation, then 
the algebraic rate is 0(d~ k )] 

• if G(t) is analytic, then there exists a positive p such that the algebraic rate is 

o( P - d ). 



This imposes some minimal requirements for A(t) and B(t) for our numerical approach to 
work. Indeed, assuming A(t) of bounded variation and B(t) differentiable with derivative 
of bounded variation guarantees that the error converges to zero for large approximation 
orders. 

We now show that if e d — » as d — > oo, we can build a hierarchy of moment relaxations of 
Q involving only approximate data and converging to p*. For this, let F d {t) denote the 
Chebyshev interpolant of F(t) on d points, and let h d := F d (t F )x(t F ) — F~[ 1 (t I )x(ti). 



Lemma 3 Consider the following relaxed moment problem with approximate data: 

p* d = min \\p\\tv 

s.t. \$G%dfi-h d \ < e d (\\B\\ \x(t F )\ + \\p\\ T v)- 

Then p* d t P* as d ^ oo 



Proof: By the same argument as Lemma [TJ a solution for (J8J) is attained. Furthermore, 
any solution p* of Q is feasible for ([8]), as 

\jG T d dp*-h d \ < \jG T dp* -h\ + \j(G T d -G T )dp*\ + \\B\\\x(t F )\\h d -h\ 
< + e d ||/i*|| T v + erf||-B|| \x{t F )\. 

Therefore, p* e < p* holds. 

For the convergence, consider the auxiliary relaxed problem with the true data instead: 

P*d= min IHlTV /qx 

s.t. \jG T dp-h\<2e d (\\B\\\x(t F )\ + \\p\\ TV ). {J> 

By similar arguments, we have that p* d < p* d < p* . Because of this uniform bound on the 
minimizers of (|9]), standard arguments (see for instance the second part of the Banach- 
Saks-Steinhaus theorem in [191 Section 13.5]) show that p* < liminf p* d1 which show indeed 
that p* d -> p*, hence p* d -> p*. □ 
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6 Solving the LP on measures 



To summarize, we have formulated our LTV optimal control problem as a primal-dual 
LP pair ([5j(6]). The coefficients in these LP problems are entries of matrix function G(t) 
and vector h. These data are calculated by numerical integration, using Chebyshev poly- 
nomial approximations, as explained in the last section. For a given degree d of the 
polynomial approximation, LP problem ^ is a particular instance of a generalized prob- 
lem of moments. It can be seen as an extension of the approach of [11 J which was originally 
designed for classical optimal control problems with polynomial dynamics. Alternatively, 
it can also be understood as an application of the approach of [6J, but after integration 
of the ODE, which is here possible because of linearity of the dynamics in the state and 
control. 

An infinite-dimensional LP on measures can be solved approximately by a hierarchy of 
finite-dimensional linear matrix inequality (LMI) problems, see [HJ EJ [9] for details (not 
reproduced here). The main idea behind the hierarchy is to manipulate each measure via 
its moments truncated to degree d. Note that here the measures are univariate (depending 
on time only). Therefore, for a problem involving polynomials up to degree d, the d-th 
LMI condition is necessary and sufficient. The hierarchy presented in this paper comes 
thus from the polynomial approximation of the continuous data, not from the moment 
truncation. In contrast, when dealing with multivariate measures as in [UJ [BJ [§], we use 
a hierarchy of necessary LMI conditions which become sufficient only asymptotically. 

Generally speaking, the number of variables in an LMI relaxation of order d of a multi- 
variate measure LP grows polynomially in d, but the exponent is the number of variables 
entering the measures. If the number of variables is equal to 5 or more, the growth is 
fast, and only LMI relaxation of small orders can be solved at a reasonable computational 
cost. It is therefore crucial to reduce as much as possible the number of variables entering 
the measures, so as to reduce the overall computational burden. One contribution of our 
paper is precisely to show that for LTV optimal control problems, we can manipulate 
measures of the time variable only. This allows for LMI relaxations of large order to be 
solved, with d a few hundreds. 

7 Examples 

7.1 Scalar polynomial example 

Consider the optimal control problem Q: 

p* = inf ||m||i 

s.t. x(t) 
x(0) 




u(t)\dt 



= t(l -t)u(t) 
= 0, x(l) = 1 
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where the minimization is w.r.t. a function u G L 1 ([0, 1]). Since -^ 1 ([0, 1]) is not the dual 
of any normed space, the problem must be relaxed to the optimal control problem (|2]): 

p* = min ||/x||tv := / \n\{dt) 

Jo 

s.t. x(dt) = t(l - t)n(dt) 
ar(0) = 0, x{\) = 1 

where the minimization is now w.r.t. a measure \i G M([0, 1]). Following the approach 
of Section |3j we readily obtain F(t) = 1, G(t) = t(l — t) and hence the moment problem 
fl: 

p* = min ||/x||tv 

s.t. (t(l -t),fj) := / t(l-t)n(dt) = l. 

Jo 

Decomposing fi = fi + — jjT as a difference of two nonnegative measures we obtain the LP 
(§: 

p* = inf (l,n+) + (l,fi-) 

s.t. (*(1 -t),fi + ) - (t(l -t),fT) = 1 
Ai + > 0, //- > 

where the minimization is w.r.t. measures and the LP @: 

d* = sup y 

s.t. -1 < *(1 - t)y < 1, Vte[0,l] 

where the maximization is w.r.t. For this latter problem, we readily obtain that 

the maximum d* = 4 is attained for the choice y = 4. Since polynomial 4i(l — t) attains 
its maximum at t — \ it follows from the discussion after Lemma [2] that the optimal 
measure solution is fi + = 45 1 and fi~ = 0, achieving p* = 4. 

As this problem has polynomial data, we can readily find the optimal solution, with the 
following GloptiPoly 3 script |8] , instead of using the approximation techniques of Section 
V: 



>> mpol tp tn; 

>> mp = meas(tp); °/„ measure \mu~ + 
>> mn = meas(tn); °/„ measure \mu~- 
>> P = msdp (min (mass (mp) +mass (mn) ) , ... 
mom(tp*(l-tp))-mom(tn*(l-tn)) == 1, ... 
tp*(l-tp) >= ... °/ support of \mu~+ 
tn*(l-tn) >= 0); % support of \mu~- 
>> [stat, obj] = msol(P); 
obj = 

4.0000 

>> double (mmat (mp) ) % moment mat. of \mu~+ 
ans = 

4.0000 2.0001 

2.0001 1.0001 
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>> double (mmat (mn) ) % moment mat. of \mu~- 
ans = 

1.0e-08 * 

0.2673 0.2036 

0.2036 0.2186 

7.2 Nonsmooth trajectories 

This example is meant to illustrate the numerical difficulty we may face when integrating 
LTV ODEs with discontinuous state matrix A(t). Consider the following optimal control 
problem Q: 

p* = inf \\u\\i 

s.t. x(t) = sign(t)x(t) + u{t) 
x\-l) = -1, x(l) = 1 

where the minimization is w.r.t. a function u G 1,1]). Upon solving the ODE 

F{t) = sign(t)F(t), F(-l) = 1 we obtain F(t) = e" 1+|t| , G(t) = e 1 " 1 * 1 and hence the 
moment problem Q: 

p* = inf ||jt/||rv 

s.t. (e 1 -l*l,^) = 2 

where the minimization is now w.r.t. a signed measure \i G M ([— 1, 1]). The dual problem 
^ reads: 

d* = sup 2y 

s.t. — 1 < 2/e 1 1*1 < 1 VtG[-l,l] 

where the maximization is w.r.t. a real scalar y. The optimal dual solution can easily be 
found analytically as y* = e" 1 , from which it follows that the optimal primal solution is 
li* = 2e- 1 5 t=0 . 



d 


\\G — Gd\ oo 


\P* - Pd\/\p*\ 


10 


0.1995 


0.1243 


20 


0.0899 


0.0555 


30 


0.0579 


0.0357 


40 


0.0427 


0.0263 


50 


0.0338 


0.0208 


60 


0.0280 


0.0172 


70 


0.0239 


0.0147 


80 


0.0208 


0.0128 


90 


0.0184 


0.0114 


100 


0.0166 


0.0102 



Table 1: Approximation errors vs. polynomial approximation orders for discontinuous 
dynamics. 

If we want to apply the numerical integration approach of Section [5j discontinuity of 
A(t) = sign(t) turns out to be a problem. Indeed, we can check that G(t) presents a cusp 
at t — 0. That is, for an even number d of Chebyshev interpolation points, the maximum 
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A(t): 



B(t) = 



n e cos v 



it) 



In e sin u{t J 






' 1 + e cos )^(t) \ 3 
, 1 - e 2 / 
/ 1 + e cos \ 3 






" 








_2 
n 








n 2 



-2 n e sin 



' 1 + e cos i^(t) 1 



n (3 + e cos v(t)) 



' 1 + e cos f(t) \ 3 



_ (1 + e cos f(t)r 
-2n 

h - r 2 ^ 3 / 2 



e cos f 



(l- e 2 ) 3 / 2 




(10) 



interpolation error ||G — G*d||oo is located where the optimal impulse should be, while for 
an odd d the interpolant is exact at the cusp. In Table [T] we present the approximation 
errors on function G and on the optimal cost p* d as functions of d, found by our numerical 
method. The results exhibit the expected linear decrease in the approximation error, 
resulting in a linear decrease in the optimal cost. The error for moderate orders (d ~ 100) 
is acceptable though far from excellent. 



d 



2 

4 

6 

8 

10 

12 

14 

16 

18 



\G — G n 



2.1187- lO" 1 
1.0851 • 10~ 3 
2.2540 ■ 
2.5115 
1.7422 • 
1.1336 • 
< 10 



■ 10" 6 

■ io- 9 
io- 12 
io- 15 

-15 

< 10~ 15 

< 10~ 15 



\p* -Pd\/\p* 



3.4307 
1.0406 
1.5281 
1.3166 
5.9438 
5.8093 
1.3729 
3.1846 
2.5462 



• IO" 3 
■ 10~ 5 
•IO" 8 

• io- 11 

• io- 13 

• io- 12 

• io- 11 

• io- 12 

• io- 10 



Table 2: Approximation errors vs. polynomial approximation orders for piecewise analytic 
dynamics. 

Now, if we split the time interval around the cusp, the restrictions of G(t) on the intervals 
[— 1, 0] and [0, 1] are respectively Gi{t) := e 1+t and Gn(t) := e 1_ * which are both analytic. 
The updated moment problem with piecewise analytic data reads: 



P 



inf 



s.t. 







\\tv + ||/4r||tv 
G L (t)ii L {dt) + 



-i 



G R (t)fj, R (dt) 



In Table [2j we present the updated polynomial approximation errors, which show the 
expected exponential decrease. For the cost, the decrease is exponential until a relative 
error of about 10 -10 after which the error comes from the semidefmite solver (used to solve 
the LMI hierarchy of the moment problem), not from the polynomial approximation. 



7.3 A fuel-optimal linear impulsive guidance rendezvous prob- 
lem for an elliptic reference orbit 

Finally, an illustration based on a realistic case of a far range rendezvous in a linearized 
gravitational field is given. The general framework of the minimum-fuel fixed-time copla- 
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nar rendezvous problem in a linear setting is recalled in [2], where an indirect method 
based on primer vector theory is proposed. Under Keplerian assumptions and for an 
elliptic reference orbit, the complete rendezvous problem may be decoupled between the 
out-of-plane rendezvous problem for which an analytical solution may be found and the 
coplanar problem. Therefore, only a coplanar elliptic rendezvous problem based on the 
Tschauner-Hempel equations [24] is studied thereafter. The associated optimal control 
problem has a 4-dimensional (n = 4) state vector composed of the relative posi- 
tions (denoted here x(t), z(t)) and respective velocities in the LVLH frame [H] and a 
2-dimensional (m = 2) control vector u(t) (one control in the x-direction, one control in 



the ^-direction). The state space matrices A{t) and B{t) are given in equ. (10), with 
n = 34.0094 is the mean angular motion, e = 4.0000 ■ 10 -3 is the eccentricity, u(t) is the 
true anomaly of the reference orbit satisfying for all t > the Kepler equation: 

nt = E(t) - e sinE(t) 

v{t) A + eV /2 E(t) 
tan = tan — ^ 



where E(t) is the eccentric anomaly. 

For this problem, G(t) can be approximated below the 10 -8 resolution of the SDP solver 
by polynomials of degree 100 on the given time interval. This implies that the problem 
can be solved numerically by LMI relaxations of order 50, with a computational load of 
a few seconds. The assembly of G(t) using Chebfun routines take a few seconds as well. 
All of this was done without any sort of problem-specific optimization, and with standard 
Matlab code. 

To solve the optimal control problem, direct methods based on the solution of a linear 
programming (LP) problem can be used as in [161 H3]- F° r an a priori fixed number of 
impulsive maneuvers at given times, an LP problem is formulated and solved numerically. 
Its solution is therefore suboptimal, depending strongly upon the number of impulsions. 
This makes it hard to evaluate how far the solution could be from the global optimum. 

The particular instance, studied in this paper, is borrowed from the PRISMA test bed 
and GNC experiments from [4]. PRISMA programme is a cooperative effort between the 
Swedish National Space Board (SNSB), the French Centre National d'Etudes Spatiales 
(CNES), the German Deutsche Zentrum fiir Luft- und Raumfahrt (DLR) and the Dan- 
ish Danmarks Tekniske Universitet (DTU) [10]. Launched on June 15, 2010 from Yasny 
(Russia), it was intended to test in-orbit new guidance schemes (particularly autonomous 
orbit control) for formation flying and rendezvous technologies. This mission includes 
the FFIORD experiment led by CNES, which features a rendezvous maneuver (formation 
acquisition). To save fuel and allow for in-flight testing throughout the FFIORD exper- 
iment, the rendezvous maneuver must last several orbits. Duration of the rendezvous is 
approximately 14.25 orbital periods, each of duration 5920s, i.e. tj = and tF = 84360 
in problem Q. 

The moment LMI approach is compared to the classical 2-impulse solution and to the 
direct approach with 20 and 2000 pre-assigned evenly distributed impulses. These re- 
sults are presented in Table [3j where tk are the impulsion times and w(^fc) the impulsion 
directions, for k = 0, 1, 2, 3. It is interesting to note that the 2-impulse solution has a 
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2-impulse 


LP (20 i.) 


LP (2000 i.) 


T 1\ IT 

LMI 


to 














h 




8880 


2321 


2140 


t 2 


- 


75480 


2363 


82350 


t 3 


84360 


84360 


82334 




u(t ) 


0.02 


0.0006 


0.0174 


0.0172 


-0.1985 















0.0014 



r\ A 1 AO 

0.0193 



0.0011 



u(t 2 ) 




-0.0146 



-0.0003 



-0.019 



u(t 3 ) 


-0.0215 
0.218 


-0.0207 



-0.019 





cost m/s 


0.4588 


0.04072 


0.03736 


0.03736 



Table 3: Comparisons of results of moment LMI approach, two-impulse solution and 
direct LP solution for 200 and 2000 impulses for the PRISMA case study 

prohibitive cost when compared to the optimal solution found by the direct method and 
the moment LMI approach. This is mainly due to the extra thrusts in the ^-direction 
that are not necessary to realize this particular rendezvous. Indeed, a remarkable feature 
of the optimal solution is that it exhibits a terminal coast, unlike the 2-impulse and di- 
rect 20-impulse solution. Note also that the 2-impulse solution leads to a very different 
trajectory in the orbital plane as shown by Figure [TJ The optimal fuel-cost computed via 
the moment LMI approach is less than half the one obtained in [I] via a pseudo-closed 
loop technique (0.086 m/s). Figure [2] shows the importance of the pre-assigned number 
of impulses for the direct method to recover the optimal solution. Indeed, for a small 
number of impulses, the solution given by the direct method is a crude approximation of 
the optimal solution obtained by the moment LMI approach and by the direct approach 
with 2000 pre-assigned impulses. The impulse locations are indicated on the trajectories 
(red triangles for the optimal solution and red squares for the 20-impulse solution). The 
differences of locations of the impulses for the four methods used on the PRISMA example 
are depicted on Figure [3j 

8 Conclusion 

In this paper, we revisit classic impulsive control theory from the 1960s with modern 
numerical tools stemming from convex programming and approximation theory. The 
theory leads to the reformulation of a control problem as a one- dimensional problem of 
moments. The numerical tools allow for its efficient numerical solution without any expert 
knowledge needed. 

Several extensions of our work are possible and will be included in an extended version of 
this paper. First of all, the class of dynamics can easily be extended to cover those with 
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Figure 1: Trajectories in the orbital plane: 2-impulse solution (pink), optimal solution 
(blue), Direct 20- impulses solution (cyan). 



a forced drift term, i.e. to dynamics 



x = A(t)x{t) + B(t)u(t) + w(t) 

with w(t) e C 1 ([t/, tp], M. n ). The second main extension is to enforce positivity constraints 
on linear combinations of the states using additional positive measures. Finally, the 
method can also consider other norms of the form ||w(-)|| p := f\u(t)\ p dt with p = 2,3, 
instead of the p — 1 case that is the focus of this paper. This can be done since norms 
with integer exponents are semidefmite representable. 



For the specific case of orbital rendezvous as exemplified in Section 7.3, one can exploit 



the specific symmetries of (10) to approximate G(t) off-line on half an orbital period, 
and partition the time interval accordingly such as presented in Ex. |7.2| This will lead 



to accurate, low order polynomial approximations, and one could expect to obtain com- 
putational loads compatible for an online implementation of a model predictive control 
loop. 
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Figure 2: Trajectories in the orbital plane: optimal solution (blue), Direct 20-impulses 
solution (cyan) . 
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